Singular Value Decomposition#
Background#
Many undergraduate physics majors take at least one mathematical methods course, in which the curriculum usually includes material on the important linear algebra concept of eigensystems; namely, the determination of eigenvectors and eigenvalues.
Eigensystems appear in all disciplines of physics, from classical mechanics to quantum mechanics; thus, while it is generally expected that a physicst can determine a system’s eigenvectors and eigenvalues, practical application is sometimes unintuitive.
This module focuses on the singular value decomposition (SVD), which decomposes an \(m \times n\) rectangular matrix \(A\) into constitutent matrices,
\begin{equation} A = U S , V^{T} \end{equation}
where \(U \in \mathbb{R}^{m \times m}\) and \(V \in \mathbb{R}^{n \times n}\) are orthogonal, and the matrix \(S \in \mathbb{R}^{m \times n}\) is a diagonal matrix of the singular values of \(A\).
There exist some differences between a canonical eigendecomposition and the SVD, wherein eigendecomposition via diagonalization of a square matrix \(A\) is of the form
\begin{equation} A = W 𝛬 , W^{-1} \end{equation}
in which the columns of \(W\) are the eigenvectors of \(A\), and the diagonal matrix \(𝛬\) holds the corresponding eigenvalues. Note that while eigendecomposition only works for square matrices, SVD is computable for rectangular matrices.
The relationship between SVD and diagonalization is that the SVD allows one to determine the eigenvectors and eigenvalues of \(AA^{T}\) (for which the columns of \(U\) are the eigenvectors) or \(A^{T}A\) (for which the columns of \(V\) are the eigenvectors); the matrix \(S\) thus corresponds to the square roots of the non-zero eigenvalues of either permutation of the product of \(A\) and its transpose.
Ironically, the problems that SVD is adept at solving are not readily apparent; however, this algorithm is crucial for a litany of applications, such as:
Image compression
Principal component analysis (PCA)
Clustering
Dimensionality reduction
Linear system solvers
This list is by no means exhaustive, but it does serve to emphasize the power of the SVD with respect to manipulating data. As is the case for many linear algebra routines, the SVD algorithm is computationally expensive when poorly implemented and fairly complicated when implemented efficiently. Fortunately, the importance of SVD in data science workflows has seen it implemented in a multitude of packages across many languages - for this module, Python’s numpy/scipy implementions will be utilized.
This module will attempt to explore applications of SVD in image compression, clustering, and PCA.
#@title Solutions
import os
import hmac
import hashlib
def encrypt_solution(solution):
"""
Function for unidirectional encryption of solution strings
"""
salt = os.urandom(16)
hash = hashlib.pbkdf2_hmac('sha256', solution.encode(), salt, 10000)
return (salt, hash)
def check_answer(solution, salt, hashed_guess):
check = hmac.compare_digest(
hashed_guess, hashlib.pbkdf2_hmac('sha256', solution.encode(), salt, 10000)
)
return check
possible_solutions = {
'exercise13' :
[
(b'\x1e\x8a\x012|\x88x\xb0Q\x05\x87\x9b6\x95?\xd5', b'\x80\xd8&\xbc7\xd9d\xd5~x\xf4W6\x06h5\xfe\tE\x11T\xa9\xa97i6D#\xb1K\xb6\xb8'),
(b'\x95W\xe5\xf0\x0e\xfc<\xf9\x03\xdc\xdc\x9e\x16\xea\xe7\xca', b'&\x16\xc7]\xaa\xb5\xbb\x8aY\x82\xe4\x86\nL\x81\x07cG\xcb\xfe\x96\x05\xe0\x19%;Sf\xc6\xb2\xf9\x87'),
(b'\xdc\xca\xa8\xb74\x8e\xf4\x0c\xd5\xe2\x9eHJ`\xce\xe4', b'A5\xbaSu\xb6M\xcbT\xfa\xcc\x86\xa1`\x18\x0c_\xef[\xbf\x03;\xd68\xd3-\xd8\xa1\xb3\xf2q\xbc'),
(b'SB\xf9\x81Y\x89\xc9iF\xadI?\xf5\x91V\xb3', b'X5\xc6M\x9d\xda\xbc!A\x9f\x1d\xd2\xac\x9dH[\xe5k\xd3\x14\xc5\xa2c`\xde\xa4M\xf36\x06kE'),
(b'#\xea\xe3f\xe6\xd9N\x92m\x1cH\xec\x87\xcb\xc3N', b'\x84\x9e\xf9]\xc2\xd7&\x18L\xa7\xdc\xad\x8dh\xcf\x7f\x82H\x96\xde\x04\xe3\xa62\xcbzP\x0f\xb4\x92We'),
(b'Z\x15t\x18n:VO\x19-\xbc\x93\xce\xf2g\xc4', b'b\xad\xe3<\xe8s\xa5\x15\xb8\x94\xdfX?\xdf\x07\r\xde\xaeB\x95\xe6^k\x89\xfb(\xadJ\x8fcA2'),
(b'\xd9P\x14\x10\x0e\xf3M\xc7\x8c\xf4\xae]\x83y\xb7\x8a', b'\xb2*l\xcb$oH\x9b\xfe>\x15!o\xd6\xf6\x83$w.\xf3ju\xaf\xb4I\x8c\x91\xa5\x01\xd1j\x9b'),
(b'\xc8\xdc\xeb\xf7\x82RW\xcf\xe3>\xd7C\xc0h\x1bb', b'\x15\x9a\xdc6z7\xdb@{\x8aA\xbb\xd5\x9e\x9eP\xe7\x93\x0c\x01T\xe8x\x7f\x81\xe2\x9e\xe1.\xaa\xf7U'),
(b'\xe3\xa0\xfe\xa9\x8d\xb6\x0cl\xa2Cl\xfbJ:\xf7d', b'\x18\xf2-\xd8%\xd9Q\xd6/r\xa8b""`\xd2VN\xa9O$\x13\xa9\xd0kX\n`\xfdd\xa5\x12')
]
}
Image Compression#
With the recent launch of the James Webb Space Telescope (JWST), cosmological images of unprecedented quality are being regularly transmitted back to Earth from 1.5 million kilometers away. With a radio capable of 28 Mbps throughput, the bandwidth of the JWST transmitter rivals and even surpasses that of some internet connections on Earth. However, despite possessing a datalink with speeds comparable to broadband internet, the telescope only has enough onboard storage for approximately 65 GB of data.
From the perspective of the technical limitations of an apparatus as advanced as the JWST, the importance of good data compression is apparent, potentially shrinking a hundred megabyte image into one of only tens of megabytes.
Many modern compression algorithms are considered lossless, such as the algorithms used by the JWST, in that they are able to reconstruct the compressed data into its original representation with no loss of fidelity. However, many lossy algorithms, such as the ubiquitious JPEG image format, are suitable nonetheless for everyday use and are still the standard for many technical applications.
This section will investigate the feasibility of using SVD to reduce the amount of data needed to render an image with acceptable quality.
#@title Application
#@markdown Suppose there exists a telescope installation on a remote island tasked solely with the purpose of capturing images of Saturn. One evening, a particularly stunning image is captured and the scientists on the island wish to send it to their colleagues.
#@markdown
#@markdown ---
#@markdown *Execute this cell to retrieve image.*
import numpy as np
import plotly.express as px
import imageio as iio
import xarray as xr
# Despite the description above, this image was captured by NASA from a telescope in space.
img = iio.imread('https://images-assets.nasa.gov/image/PIA05982/PIA05982~medium.jpg')
fig = px.imshow(img)
fig.update_layout(showlegend=False)
fig.update_xaxes(visible=False)
fig.update_yaxes(visible=False)
fig.update_layout(autosize = False, width = 1280, height = 663)
fig.show()
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In [2], line 7
1 #@title Application
2 #@markdown Suppose there exists a telescope installation on a remote island tasked solely with the purpose of capturing images of Saturn. One evening, a particularly stunning image is captured and the scientists on the island wish to send it to their colleagues.
3 #@markdown
4 #@markdown ---
5 #@markdown *Execute this cell to retrieve image.*
6 import numpy as np
----> 7 import plotly.express as px
8 import imageio as iio
9 import xarray as xr
ModuleNotFoundError: No module named 'plotly'
Unfortunately, there is no internet on the island and all communication must be transmitted via satellite; however, the satellite is only over the island for a short period of time during which uplink may occur. This limits the size of transmitted data, since any uploads must be completed during a brief window of time. Ultimately, this doesn’t necessarily make the image significantly worse and much of the lost data may be compensated for by other means, such as interpolation/upsampling to recover lost dimensionality and denoising to clean up any artifacts introduced during compression.
At its original size, the image is too large to fully transmit before the satellite is out of range. However, consider that an image is simply a matrix of pixel intensities; thus, an image is a matrix \(A\) which may be decomposed into its singular values. SVD may be used to reduce the amount of data needed to represent the image, making it possible to upload the image to the satellite.
Implementation#
Singular value decomposition is implemented in Numpy as
u, s, v = np.linalg.svd(A)
Numpy’s implementation of SVD does not natively work on RGB images, since they exist as a 3D array containing a slice for each channel of the image. The code below shows how np.linalg.svd may be extended to work on RGB images.
# `np.linalg.svd` does not perform SVD on RGB images. The following function takes
# an input `image` and a desired number of singular values `n` and computes a
# channel-wise SVD, after which the channels are collated and re-formed into
# a uint8 ndarray representable as an RGB image.
def rgb2svd(image, n):
"""
Compute the SVD of a 3-channel RGB image.
This function additionally takes a parameter `n` to specify how many
singular values should be extracted from the image, since this function takes
significantly longer to execute for large images.
...
Parameters
----------
image : ndarray
Multi-channel matrix containing image data.
n : int
Number of singular values to extract from the image.
"""
# Instantiate a list to accumulate the R,G,B channels of the image.
acc = []
# Iterate through the 3 channels of the image and decompose each one.
for i in range(image.shape[2]):
# This line performs the actual SVD on the current channel denoted by `i`.
u, s, v = np.linalg.svd(image[:,:,i], full_matrices=False)
# Recompose the image using `n` singular values from the SVD.
# Necessary to constrain pixel intensities to be positive due to uint8 cast.
aa = np.abs(np.dot(u[:,:n], np.dot(np.diag(s[:n]), v[:n,:])))
# Rescale the RGB values between [0, 255] and cast to uint8.
res = (aa * 255 / np.max(aa)).astype('uint8')
# Collect the compressed channel in `acc`.
acc.append(res)
# Create the new, compressed image by stacking the channels in the third dimension.
return np.dstack(tuple(acc))
Generally, it is beneficial to pre-compute large amounts of data rather than computing on-the-fly inside a loop. This function creates a “gallery” of images reconstructed from increasing numbers of singular values, which is then plotted below. Use the start, stop, and step parameters below to adjust how the singular values are sampled.
# Alter these parameters to change the number of singular values sampled from the image.
# start - the number of singular values at which to begin sampling
# stop - the number of singular values at which to stop sampling
# step - the singular value sampling rate
start, stop, step = 1, 120, 10
# This function may take a while to run depending on the size of the input image.
# Change the `step` parameter above to reduce the runtime of this function.
def svd_gallery(image):
"""
Generate a 'gallery' of SVD compressed images.
...
Parameters
----------
image : ndarray
Multi-channel matrix containing image data.
"""
gallery = []
x, _, _ = image.shape
for i in np.arange(start, stop, step, dtype=int):
if i == 0: i = 1
gallery.append(rgb2svd(image, n=i))
return np.asarray(gallery)
img_gallery = svd_gallery(img)
#@title Visualization of singular value image reconstruction
#@markdown The plot below contains a slider object, which may be used to cycle through the reconstructed images - from the image reconstructed from the fewest singular values to an image reconstructed from significantly more. <br>
#@markdown The number of singular values is truncated due to the size of the data.
#@markdown Use the slider in the plot below to visualize the reconstructed images and determine at what point the diminishing returns on the number of singular values used to reconstruct the image become apparent.
# Plotly Express requires `xarray` data to generate animations.
h, w, d = img.shape
sval_coord = np.fromiter(
(1 if i == 0 else i for i in range(start, stop, step)),
dtype=int)
y_coord = np.arange(h)
x_coord = np.arange(w)
d_coord = np.arange(d)
dsimg = xr.DataArray(
data = img_gallery,
dims = ['Number of singular values sampled by reconstructed image', 'y', 'x', 'channel'],
coords = [sval_coord, y_coord, x_coord, d_coord],
attrs= {
'long_name': 'Reconstructed Images From N Singular Values',
}
)
fig2 = px.imshow(
dsimg,
animation_frame="Number of singular values sampled by reconstructed image",
binary_string=True,
title=dsimg.attrs["long_name"])
#fig2["layout"].pop("updatemenus")
fig2.update_layout(showlegend=True)
fig2.update_xaxes(visible=False)
fig2.update_yaxes(visible=False)
fig2.update_layout(autosize = False, width = 1280, height = 663)
fig2.show()